Load the required packages and set a seed for the random number generator.

library(tidyverse)

library(rstan)
# set cores to use to the total number of cores (minimally 4)
options(mc.cores = max(parallel::detectCores(), 4))
# save a compiled version of the Stan model file
rstan_options(auto_write = TRUE)

library(brms)

# install faintr with devtools::install_github('michael-franke/bayes_mixed_regression_tutorial/faintr', build_vignettes = TRUE)
library(faintr)

set.seed(123)

Read the data:

raw_data = read_csv("results.csv")
## Warning: 1470 parsing failures.
##  row      col           expected                                                                                                                                                                                                                      actual          file
## 7106 comments 1/0/T/F/TRUE/FALSE I only assigned values between 2 and 5 because I expected nicer and uglier pictures to come. My values are not correct as absolute ratings. They are influenced by relative comarision to what appeared immediately before. 'results.csv'
## 7107 comments 1/0/T/F/TRUE/FALSE I only assigned values between 2 and 5 because I expected nicer and uglier pictures to come. My values are not correct as absolute ratings. They are influenced by relative comarision to what appeared immediately before. 'results.csv'
## 7108 comments 1/0/T/F/TRUE/FALSE I only assigned values between 2 and 5 because I expected nicer and uglier pictures to come. My values are not correct as absolute ratings. They are influenced by relative comarision to what appeared immediately before. 'results.csv'
## 7109 comments 1/0/T/F/TRUE/FALSE I only assigned values between 2 and 5 because I expected nicer and uglier pictures to come. My values are not correct as absolute ratings. They are influenced by relative comarision to what appeared immediately before. 'results.csv'
## 7110 comments 1/0/T/F/TRUE/FALSE I only assigned values between 2 and 5 because I expected nicer and uglier pictures to come. My values are not correct as absolute ratings. They are influenced by relative comarision to what appeared immediately before. 'results.csv'
## .... ........ .................. ........................................................................................................................................................................................................................... .............
## See problems(...) for more details.
glimpse(raw_data)
## Observations: 27,685
## Variables: 29
## $ submission_id <dbl> 605, 605, 605, 605, 605, 605, 605, 605, 605, 605, …
## $ QUD           <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ RT            <dbl> 15344, 8314, 65189, 10070, 5679, 2714, 4583, 2000,…
## $ age           <dbl> 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23…
## $ artist        <chr> NA, NA, NA, "P", "B", "G", "G", "P", "P", "P", "G"…
## $ color         <chr> NA, NA, NA, "colored", "monochrome", "colored", "c…
## $ color2        <chr> NA, NA, NA, "sepia", "monochrome", "colored", "col…
## $ comments      <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ date          <dbl> NA, NA, NA, 1910, 1911, 1915, 1912, 1911, 1911, 19…
## $ education     <chr> "Graduated College", "Graduated College", "Graduat…
## $ endTime       <dbl> 1.563018e+12, 1.563018e+12, 1.563018e+12, 1.563018…
## $ expected      <chr> "74", "7", "FELOPZD", NA, NA, NA, NA, NA, NA, NA, …
## $ experiment_id <dbl> 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77…
## $ gender        <chr> "female", "female", "female", "female", "female", …
## $ languages     <chr> "German", "German", "German", "German", "German", …
## $ min_chars     <dbl> 0, 0, 1, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ option1       <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ option2       <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ optionLeft    <chr> NA, NA, NA, "not at all", "not at all", "not at al…
## $ optionRight   <chr> NA, NA, NA, "very", "very", "very", "very", "very"…
## $ picture       <chr> "vision/ishihara.png", "vision/ishihara2.jpg", "vi…
## $ pictureNumber <dbl> 0, 0, 0, 32, 65, 75, 84, 56, 71, 120, 89, 46, 16, …
## $ question      <chr> "Which number do you see here?", "Which number do …
## $ response      <chr> "74", "7", "FELOPZD", "4", "5", "3", "4", "3", "3"…
## $ startDate     <chr> "Sat Jul 13 2019 13:21:37 GMT+0200 (CEST)", "Sat J…
## $ startTime     <dbl> 1.563017e+12, 1.563017e+12, 1.563017e+12, 1.563017…
## $ timeSpent     <dbl> 13.91613, 13.91613, 13.91613, 13.91613, 13.91613, …
## $ trial_name    <chr> "colourVisionTest1", "colourVisionTest2", "visionT…
## $ trial_number  <dbl> 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13…

Rename variables

data = rename(raw_data,
  submissionID = `submission_id`,
  trialName = `trial_name`
)

glimpse(data)
## Observations: 27,685
## Variables: 29
## $ submissionID  <dbl> 605, 605, 605, 605, 605, 605, 605, 605, 605, 605, …
## $ QUD           <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ RT            <dbl> 15344, 8314, 65189, 10070, 5679, 2714, 4583, 2000,…
## $ age           <dbl> 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23…
## $ artist        <chr> NA, NA, NA, "P", "B", "G", "G", "P", "P", "P", "G"…
## $ color         <chr> NA, NA, NA, "colored", "monochrome", "colored", "c…
## $ color2        <chr> NA, NA, NA, "sepia", "monochrome", "colored", "col…
## $ comments      <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ date          <dbl> NA, NA, NA, 1910, 1911, 1915, 1912, 1911, 1911, 19…
## $ education     <chr> "Graduated College", "Graduated College", "Graduat…
## $ endTime       <dbl> 1.563018e+12, 1.563018e+12, 1.563018e+12, 1.563018…
## $ expected      <chr> "74", "7", "FELOPZD", NA, NA, NA, NA, NA, NA, NA, …
## $ experiment_id <dbl> 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77, 77…
## $ gender        <chr> "female", "female", "female", "female", "female", …
## $ languages     <chr> "German", "German", "German", "German", "German", …
## $ min_chars     <dbl> 0, 0, 1, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ option1       <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ option2       <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ optionLeft    <chr> NA, NA, NA, "not at all", "not at all", "not at al…
## $ optionRight   <chr> NA, NA, NA, "very", "very", "very", "very", "very"…
## $ picture       <chr> "vision/ishihara.png", "vision/ishihara2.jpg", "vi…
## $ pictureNumber <dbl> 0, 0, 0, 32, 65, 75, 84, 56, 71, 120, 89, 46, 16, …
## $ question      <chr> "Which number do you see here?", "Which number do …
## $ response      <chr> "74", "7", "FELOPZD", "4", "5", "3", "4", "3", "3"…
## $ startDate     <chr> "Sat Jul 13 2019 13:21:37 GMT+0200 (CEST)", "Sat J…
## $ startTime     <dbl> 1.563017e+12, 1.563017e+12, 1.563017e+12, 1.563017…
## $ timeSpent     <dbl> 13.91613, 13.91613, 13.91613, 13.91613, 13.91613, …
## $ trialName     <chr> "colourVisionTest1", "colourVisionTest2", "visionT…
## $ trial_number  <dbl> 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13…

Select only the variables that are important for the analysis:

data = select(data, submissionID, trialName, 
         response, pictureNumber, artist, date, 
         age, gender, timeSpent, color, color2)

glimpse(data)
## Observations: 27,685
## Variables: 11
## $ submissionID  <dbl> 605, 605, 605, 605, 605, 605, 605, 605, 605, 605, …
## $ trialName     <chr> "colourVisionTest1", "colourVisionTest2", "visionT…
## $ response      <chr> "74", "7", "FELOPZD", "4", "5", "3", "4", "3", "3"…
## $ pictureNumber <dbl> 0, 0, 0, 32, 65, 75, 84, 56, 71, 120, 89, 46, 16, …
## $ artist        <chr> NA, NA, NA, "P", "B", "G", "G", "P", "P", "P", "G"…
## $ date          <dbl> NA, NA, NA, 1910, 1911, 1915, 1912, 1911, 1911, 19…
## $ age           <dbl> 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23…
## $ gender        <chr> "female", "female", "female", "female", "female", …
## $ timeSpent     <dbl> 13.91613, 13.91613, 13.91613, 13.91613, 13.91613, …
## $ color         <chr> NA, NA, NA, "colored", "monochrome", "colored", "c…
## $ color2        <chr> NA, NA, NA, "sepia", "monochrome", "colored", "col…

Exclude participants from the data that we don’t want to have. Also find out why we had to exclude them.

data_ex <- data

fehlerCV1 = 0
fehlerCV2 = 0
fehlerV = 0
fehlerExp = 0
fehlerT = 0
fehlerXP = 0


for(i in 1:(nrow(data) / 245)) {
  
  # colorVisionTest1
  if(data[((i - 1)*245 + 1), ]$response != "74") {
    data_ex = filter(data_ex, submissionID != data[((i - 1)*245 + 1), ]$submissionID)
    
    fehlerCV1 = fehlerCV1 + 1
    
  # colorVisionTest2
  } else if(data[((i - 1)*245 + 2), ]$response != "7") {
    data_ex = filter(data_ex, submissionID != data[((i - 1)*245 + 2), ]$submissionID)
    
     fehlerCV2 = fehlerCV2 + 1
     
  # visionTest
  } else if(data[((i - 1)*245 + 3), ]$response != "FELOPZD") {
    data_ex = filter(data_ex, submissionID != data[((i - 1)*245 + 3), ]$submissionID)
    
     fehlerV = fehlerV + 1
     
  # expertise 
  } else if(as.integer(data[((i - 1)*245 + 244), ]$response) > 5) {
     data_ex = filter(data_ex, submissionID != data[((i - 1)*245 + 244), ]$submissionID)
   
      fehlerExp = fehlerExp + 1
      
  # time spent 
  } else if(data[((i - 1)*245 + 1), ]$timeSpent < 5) {
     data_ex = filter(data_ex, submissionID != data[((i - 1)*245 + 1), ]$submissionID)
   
      fehlerT = fehlerT + 1
      
  # exclude participants that took part in the course Experimental Psychology Lab
  } else if(data[((i - 1)*245 + 245), ]$response == "Yes") {
    data_ex = filter(data_ex, submissionID != data[((i - 1)*245 + 245), ]$submissionID)
    
     fehlerXP = fehlerXP + 1
     
  }
  
}

glimpse(data_ex)
## Observations: 15,680
## Variables: 11
## $ submissionID  <dbl> 605, 605, 605, 605, 605, 605, 605, 605, 605, 605, …
## $ trialName     <chr> "colourVisionTest1", "colourVisionTest2", "visionT…
## $ response      <chr> "74", "7", "FELOPZD", "4", "5", "3", "4", "3", "3"…
## $ pictureNumber <dbl> 0, 0, 0, 32, 65, 75, 84, 56, 71, 120, 89, 46, 16, …
## $ artist        <chr> NA, NA, NA, "P", "B", "G", "G", "P", "P", "P", "G"…
## $ date          <dbl> NA, NA, NA, 1910, 1911, 1915, 1912, 1911, 1911, 19…
## $ age           <dbl> 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23…
## $ gender        <chr> "female", "female", "female", "female", "female", …
## $ timeSpent     <dbl> 13.91613, 13.91613, 13.91613, 13.91613, 13.91613, …
## $ color         <chr> NA, NA, NA, "colored", "monochrome", "colored", "c…
## $ color2        <chr> NA, NA, NA, "sepia", "monochrome", "colored", "col…
print(c("Number of excluded participants because of colorVisionTest1: ", fehlerCV1))
## [1] "Number of excluded participants because of colorVisionTest1: "
## [2] "18"
print(c("Number of excluded participants because of colorVisionTest2: ", fehlerCV2))
## [1] "Number of excluded participants because of colorVisionTest2: "
## [2] "2"
print(c("Number of excluded participants because of VisionTest: ", fehlerV))
## [1] "Number of excluded participants because of VisionTest: "
## [2] "24"
print(c("Number of excluded participants because of their expertise level in cubism: ", fehlerExp))
## [1] "Number of excluded participants because of their expertise level in cubism: "
## [2] "0"
print(c("Number of excluded participants because of TimeSpent: ", fehlerT))
## [1] "Number of excluded participants because of TimeSpent: "
## [2] "0"
print(c("Number of excluded participants because of participation in XPLab: ", fehlerXP))
## [1] "Number of excluded participants because of participation in XPLab: "
## [2] "5"

Calculate age mean, age range and gender distribution, number of participants:

data_ex %>% summarise(meanAge = mean(na.exclude(age)))
## # A tibble: 1 x 1
##   meanAge
##     <dbl>
## 1    29.3
gender_distribution = count(data_ex, gender) 
gender_distribution$n =  gender_distribution$n / 245

print(c("Age range: ", range(na.exclude(data_ex$age))))
## [1] "Age range: " "11"          "81"
participants = nrow(data) / 245
participants_valid = nrow(data_ex) / 245
    
print(gender_distribution)
## # A tibble: 2 x 2
##   gender     n
##   <chr>  <dbl>
## 1 female    40
## 2 male      24
print(c("Total number of participants: ", participants))
## [1] "Total number of participants: " "113"
print(c("Number of valid participants: ", participants_valid))
## [1] "Number of valid participants: " "64"

transform the data into right form:

d = filter(data_ex, trialName %in% c("ratingScaleLike", 
                                 "ratingScaleDetect")) %>%
  select(submissionID, trialName, pictureNumber, artist, date, response, color, color2) %>% 
  spread(key = trialName, value = response) 

d$ratingScaleLike <- as.numeric(d$ratingScaleLike)
d$ratingScaleDetect <- as.numeric(d$ratingScaleDetect)

d = rename(d, detectability = ratingScaleDetect, liking = ratingScaleLike)

glimpse(d)
## Observations: 7,680
## Variables: 8
## $ submissionID  <dbl> 551, 551, 551, 551, 551, 551, 551, 551, 551, 551, …
## $ pictureNumber <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,…
## $ artist        <chr> "B", "B", "B", "B", "B", "B", "B", "B", "B", "P", …
## $ date          <dbl> 1909, 1909, 1910, 1910, 1910, 1910, 1913, 1909, 19…
## $ color         <chr> "colored", "colored", "colored", "colored", "color…
## $ color2        <chr> "colored", "sepia", "sepia", "sepia", "colored", "…
## $ detectability <dbl> 6, 5, 4, 2, 3, 2, 2, 6, 2, 2, 5, 2, 4, 4, 7, 2, 1,…
## $ liking        <dbl> 2, 2, 2, 2, 2, 5, 2, 2, 2, 2, 1, 2, 2, 2, 3, 3, 4,…

Check that everything is there and of the correct type:

d_tibble <- as_tibble(d)

print(d_tibble)
## # A tibble: 7,680 x 8
##    submissionID pictureNumber artist  date color color2 detectability
##           <dbl>         <dbl> <chr>  <dbl> <chr> <chr>          <dbl>
##  1          551             1 B       1909 colo… color…             6
##  2          551             2 B       1909 colo… sepia              5
##  3          551             3 B       1910 colo… sepia              4
##  4          551             4 B       1910 colo… sepia              2
##  5          551             5 B       1910 colo… color…             3
##  6          551             6 B       1910 colo… sepia              2
##  7          551             7 B       1913 colo… color…             2
##  8          551             8 B       1909 colo… color…             6
##  9          551             9 B       1911 colo… sepia              2
## 10          551            10 P       1912 colo… color…             2
## # … with 7,670 more rows, and 1 more variable: liking <dbl>

inspect the data with x = detectability and y = liking for the different participants This part is just for our interest.

ggplot(d, aes(x = detectability, y = liking, group = submissionID)) + geom_jitter() + geom_smooth(method = "lm") + labs(x = "Detectability", y = "Liking")

inspect the data x = detectability and y = liking

ggplot(d, aes(x = detectability, y = liking)) + geom_jitter() + geom_smooth(method = "lm") + labs(x = "Detectability", y = "Liking")

ggplot(d, aes(x = detectability, y = liking)) + geom_point() + geom_smooth(method = "lm") + labs(x = "Detectability", y = "Liking")

Boxplots that shows liking and detectability

ggplot(d, aes(x = detectability, y = liking, group = detectability)) + geom_jitter() + geom_boxplot() + labs(x = "Detectability", y = "Liking")

ggplot(d, aes(x = detectability, y = liking, group = detectability)) + geom_boxplot() + labs(x = "Detectability", y = "Liking")

calculate means for detectability:

d %>% group_by(liking) %>% summarise(meanDetect = mean(detectability))
## # A tibble: 7 x 2
##   liking meanDetect
##    <dbl>      <dbl>
## 1      1       2.68
## 2      2       2.88
## 3      3       3.37
## 4      4       3.74
## 5      5       4.12
## 6      6       4.28
## 7      7       5.13

calculate means for liking:

d %>% group_by(detectability) %>% summarise(meanLike = mean(liking))
## # A tibble: 7 x 2
##   detectability meanLike
##           <dbl>    <dbl>
## 1             1     2.27
## 2             2     2.75
## 3             3     2.95
## 4             4     3.23
## 5             5     3.40
## 6             6     3.51
## 7             7     3.49

calculate means for different artists:

d %>% group_by(artist) %>% summarise(meanLike = mean(liking), meanDetect = mean(detectability))
## # A tibble: 3 x 3
##   artist meanLike meanDetect
##   <chr>     <dbl>      <dbl>
## 1 B          2.71       2.69
## 2 G          3.36       4.78
## 3 P          2.74       2.53

Boxplots and graphs that show liking and detectability (for the different artists)

dataB <- d %>% filter(artist == "B")
dataP <- d %>% filter(artist == "P")
dataG <- d %>% filter(artist == "G")

ggplot(dataB, aes(x = detectability, y = liking, group = detectability)) + geom_boxplot() + labs(x = "Detectability", y = "Liking", title = "Braque")

ggplot(dataP, aes(x = detectability, y = liking, group = detectability)) + geom_boxplot() + labs(x = "Detectability", y = "Liking", title = "Picasso")

ggplot(dataG, aes(x = detectability, y = liking, group = detectability)) + geom_boxplot() + labs(x = "Detectability", y = "Liking", title = "Gris")

ggplot(dataB, aes(x = detectability, y = liking)) + geom_point() + geom_smooth(method = "lm") + labs(x = "Detectability", y = "Liking", title = "Braque")

ggplot(dataP, aes(x = detectability, y = liking)) + geom_point() + geom_smooth(method = "lm") + labs(x = "Detectability", y = "Liking", title = "Picasso")

ggplot(dataG, aes(x = detectability, y = liking)) + geom_point() + geom_smooth(method = "lm") + labs(x = "Detectability", y = "Liking", title = "Gris")

Graphs / boxplots that show liking and detectability (for the different artists in one graph / boxplot)

ggplot(d, aes(x = detectability, y = liking, fill = artist)) + geom_point() + geom_smooth(method = "lm") + labs(x = "Detectability", y = "Liking", title = "Different artists", fill = "Artist")

ggplot(d, aes(x = as.factor(detectability), y = liking, fill = artist)) + geom_boxplot()  + labs(x = "Detectability", y = "Liking", title = "Different artists", fill = "Artist")

Graphs / boxplots that show liking and detectability (for the color vs. monochrome in one graph / boxplot)

ggplot(d, aes(x = detectability, y = liking, fill = color)) + geom_point() + geom_smooth(method = "lm" ) + labs(x = "Detectability", y = "Liking", title = "Different colors", fill = "Color")

ggplot(d, aes(x = as.factor(detectability), y = liking, fill = color)) + geom_boxplot() + labs(x = "Detectability", y = "Liking", title = "Different colors", fill = "Color")

Graphs / boxplots that show liking and detectability (for the color vs. monochrome vs. sepia in one graph / boxplot)

ggplot(d, aes(x = detectability, y = liking, fill = color2)) + geom_point() + geom_smooth(method = "lm") + labs(x = "Detectability", y = "Liking", title = "Different colors", fill = "Color")

ggplot(d, aes(x = as.factor(detectability), y = liking, fill = color2)) + geom_boxplot() + labs(x = "Detectability", y = "Liking", title = "Different colors", fill = "Color")

calculate means for different colors:

d %>% group_by(color2) %>% summarise(meanLike = mean(liking))
## # A tibble: 3 x 2
##   color2     meanLike
##   <chr>         <dbl>
## 1 colored        3.31
## 2 monochrome     2.77
## 3 sepia          2.73

Boxplots and graphs that show liking and detectability (for the different time span)

data1 <- d %>% filter(date <= 1911)
data2 <- d %>% filter(date >= 1912, date <= 1915)
data3 <- d %>% filter(date >= 1916)

ggplot(data1, aes(x = detectability, y = liking, group = detectability)) + geom_boxplot() + labs(x = "Detectability", y = "Liking", title = "Paintings from 1908-1911")

ggplot(data2, aes(x = detectability, y = liking, group = detectability)) + geom_boxplot() + labs(x = "Detectability", y = "Liking", title = "Paintings from 1912-1915")

ggplot(data3, aes(x = detectability, y = liking, group = detectability)) + geom_boxplot() + labs(x = "Detectability", y = "Liking", title = "Paintings from 1916-1919")

ggplot(data1, aes(x = detectability, y = liking)) + geom_point() + geom_smooth(method = "lm") + labs(x = "Detectability", y = "Liking", title = "Paintings from 1908-1911")

ggplot(data2, aes(x = detectability, y = liking)) + geom_point() + geom_smooth(method = "lm") + labs(x = "Detectability", y = "Liking", title = "Paintings from 1912-1915")

ggplot(data3, aes(x = detectability, y = liking)) + geom_point() + geom_smooth(method = "lm") + labs(x = "Detectability", y = "Liking", title = "Paintings from 1916-1919")

prepare the data for the analysis:

d = mutate(d, 
           liking = factor(liking, ordered = T),
           detectability = factor(detectability, ordered = T),
           color = factor(color),
           color2 = factor(color2),
           artist = factor(artist)
           )
glimpse(d)
## Observations: 7,680
## Variables: 8
## $ submissionID  <dbl> 551, 551, 551, 551, 551, 551, 551, 551, 551, 551, …
## $ pictureNumber <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,…
## $ artist        <fct> B, B, B, B, B, B, B, B, B, P, B, P, P, P, P, P, P,…
## $ date          <dbl> 1909, 1909, 1910, 1910, 1910, 1910, 1913, 1909, 19…
## $ color         <fct> colored, colored, colored, colored, colored, color…
## $ color2        <fct> colored, sepia, sepia, sepia, colored, sepia, colo…
## $ detectability <ord> 6, 5, 4, 2, 3, 2, 2, 6, 2, 2, 5, 2, 4, 4, 7, 2, 1,…
## $ liking        <ord> 2, 2, 2, 2, 2, 5, 2, 2, 2, 2, 1, 2, 2, 2, 3, 3, 4,…

ANALYSIS

Pearson correlation of liking and detectability:

cor(as.numeric(d$liking), as.numeric(d$detectability))
## [1] 0.2828698

M1. simple analysis: “objects” ratings treated interval scale / metric Calculate the difference of the sample points: compare the sample points for each category of detectability (2 to 7) to the sample points from detectability = 1.

m1 = brm(formula = liking ~ detectability,
         data = d,
         family = cumulative("logit"))

print(m1)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: liking ~ detectability 
##    Data: d (Number of observations: 7680) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept[1]       -1.75      0.03    -1.82    -1.69       3278 1.00
## Intercept[2]       -0.45      0.03    -0.50    -0.40       4259 1.00
## Intercept[3]        0.65      0.03     0.60     0.70       4905 1.00
## Intercept[4]        1.55      0.03     1.49     1.61       5275 1.00
## Intercept[5]        2.56      0.04     2.47     2.65       5246 1.00
## Intercept[6]        3.74      0.08     3.59     3.89       4742 1.00
## detectability.L     1.29      0.06     1.18     1.40       4527 1.00
## detectability.Q    -0.56      0.06    -0.67    -0.45       4243 1.00
## detectability.C     0.05      0.06    -0.07     0.16       4135 1.00
## detectabilityE4    -0.11      0.06    -0.23     0.00       4480 1.00
## detectabilityE5     0.04      0.06    -0.07     0.15       5050 1.00
## detectabilityE6    -0.05      0.06    -0.17     0.07       5236 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).
marginal_effects(m1, categorical = T)

marginal_effects(m1, categorical = F)
## Warning: Predictions are treated as continuous variables in
## 'marginal_effects' by default, which is likely invalid for ordinal
## families. Please set 'categorical' to TRUE.

M2: Model that treats detectability as monotonic and includes random effects

m2 = brm(formula = liking ~ mo(detectability) + (1|submissionID),
         data = d,
         family = "cumulative")

print(m2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: liking ~ mo(detectability) + (1 | submissionID) 
##    Data: d (Number of observations: 7680) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~submissionID (Number of levels: 64) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)     1.37      0.13     1.14     1.66        381 1.01
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept[1]       -1.12      0.18    -1.47    -0.75        107 1.01
## Intercept[2]        0.54      0.18     0.20     0.91        106 1.01
## Intercept[3]        1.89      0.18     1.55     2.26        107 1.01
## Intercept[4]        2.94      0.18     2.59     3.32        111 1.01
## Intercept[5]        4.04      0.18     3.69     4.43        113 1.01
## Intercept[6]        5.27      0.19     4.90     5.68        128 1.00
## modetectability     1.65      0.09     1.47     1.83       2145 1.00
## 
## Simplex Parameters: 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## modetectability1[1]     0.26      0.04     0.18     0.33       3212 1.00
## modetectability1[2]     0.15      0.04     0.07     0.23       2261 1.00
## modetectability1[3]     0.19      0.05     0.09     0.29       2032 1.00
## modetectability1[4]     0.17      0.06     0.06     0.28       1840 1.00
## modetectability1[5]     0.12      0.05     0.02     0.22       2353 1.00
## modetectability1[6]     0.12      0.05     0.02     0.22       2715 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).
marginal_effects(m2, categorical = T)

marginal_effects(m2, categorical = F)
## Warning: Predictions are treated as continuous variables in
## 'marginal_effects' by default, which is likely invalid for ordinal
## families. Please set 'categorical' to TRUE.

M3: Model that treats detectability as monotonic

m3 = brm(formula = liking ~ mo(detectability),
         data = d,
         family = "cumulative")

print(m3)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: liking ~ mo(detectability) 
##    Data: d (Number of observations: 7680) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept[1]       -0.66      0.05    -0.75    -0.58       3843 1.00
## Intercept[2]        0.64      0.05     0.55     0.73       4248 1.00
## Intercept[3]        1.74      0.05     1.64     1.83       4382 1.00
## Intercept[4]        2.64      0.05     2.53     2.74       4452 1.00
## Intercept[5]        3.64      0.06     3.51     3.76       4942 1.00
## Intercept[6]        4.82      0.09     4.64     4.99       5616 1.00
## modetectability     1.59      0.07     1.46     1.71       3802 1.00
## 
## Simplex Parameters: 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## modetectability1[1]     0.45      0.03     0.38     0.52       4558 1.00
## modetectability1[2]     0.17      0.04     0.09     0.26       3698 1.00
## modetectability1[3]     0.20      0.05     0.11     0.29       5008 1.00
## modetectability1[4]     0.10      0.05     0.01     0.19       5247 1.00
## modetectability1[5]     0.05      0.04     0.00     0.13       5565 1.00
## modetectability1[6]     0.03      0.02     0.00     0.08       6052 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).
marginal_effects(m3, categorical = T)

marginal_effects(m3, categorical = F)
## Warning: Predictions are treated as continuous variables in
## 'marginal_effects' by default, which is likely invalid for ordinal
## families. Please set 'categorical' to TRUE.

M4: Model that additionlly includes artist and color as an independent variable and includes random effects

m4 = brm(formula = liking ~ mo(detectability) + artist + color2 + (1|submissionID), 
         data = d,
         family = "cumulative")

print(m4)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: liking ~ mo(detectability) + artist + color2 + (1 | submissionID) 
##    Data: d (Number of observations: 7680) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~submissionID (Number of levels: 64) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)     1.41      0.14     1.16     1.71        366 1.01
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept[1]        -1.57      0.18    -1.93    -1.19        165 1.03
## Intercept[2]         0.13      0.18    -0.22     0.51        163 1.03
## Intercept[3]         1.52      0.19     1.17     1.90        165 1.03
## Intercept[4]         2.59      0.19     2.23     2.97        171 1.03
## Intercept[5]         3.71      0.19     3.35     4.10        169 1.03
## Intercept[6]         4.96      0.20     4.57     5.36        185 1.03
## artistG              0.50      0.06     0.38     0.62       3574 1.00
## artistP              0.02      0.05    -0.08     0.12       4026 1.00
## color2monochrome    -0.66      0.05    -0.76    -0.57       3390 1.00
## color2sepia         -0.32      0.07    -0.45    -0.19       3006 1.00
## modetectability      1.19      0.10     1.00     1.39       3706 1.00
## 
## Simplex Parameters: 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## modetectability1[1]     0.33      0.05     0.23     0.43       3925 1.00
## modetectability1[2]     0.15      0.06     0.04     0.27       2523 1.00
## modetectability1[3]     0.16      0.07     0.03     0.30       2284 1.00
## modetectability1[4]     0.16      0.07     0.03     0.29       3351 1.00
## modetectability1[5]     0.09      0.06     0.01     0.22       4833 1.00
## modetectability1[6]     0.11      0.06     0.01     0.24       4008 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).
marginal_effects(m4, categorical = T)

marginal_effects(m4, categorical = F)
## Warning: Predictions are treated as continuous variables in
## 'marginal_effects' by default, which is likely invalid for ordinal
## families. Please set 'categorical' to TRUE.

## Warning: Predictions are treated as continuous variables in
## 'marginal_effects' by default, which is likely invalid for ordinal
## families. Please set 'categorical' to TRUE.

## Warning: Predictions are treated as continuous variables in
## 'marginal_effects' by default, which is likely invalid for ordinal
## families. Please set 'categorical' to TRUE.

M5: Like M4, but without random effects

m5 = brm(formula = liking ~ mo(detectability) + artist + color2, 
          data = d,
          family = "cumulative")

print(m5)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: liking ~ mo(detectability) + artist + color2 
##    Data: d (Number of observations: 7680) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept[1]        -0.93      0.07    -1.07    -0.80       3577 1.00
## Intercept[2]         0.39      0.07     0.25     0.52       3776 1.00
## Intercept[3]         1.50      0.07     1.37     1.64       3919 1.00
## Intercept[4]         2.41      0.07     2.27     2.55       4035 1.00
## Intercept[5]         3.42      0.08     3.26     3.58       4208 1.00
## Intercept[6]         4.61      0.10     4.41     4.81       5103 1.00
## artistG              0.27      0.06     0.16     0.39       3489 1.00
## artistP              0.04      0.05    -0.06     0.14       5100 1.00
## color2monochrome    -0.51      0.05    -0.60    -0.42       4549 1.00
## color2sepia         -0.22      0.06    -0.34    -0.09       4492 1.00
## modetectability      1.41      0.07     1.28     1.55       3459 1.00
## 
## Simplex Parameters: 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## modetectability1[1]     0.51      0.04     0.43     0.58       4257 1.00
## modetectability1[2]     0.18      0.05     0.09     0.27       4741 1.00
## modetectability1[3]     0.18      0.05     0.07     0.28       4909 1.00
## modetectability1[4]     0.08      0.04     0.01     0.17       4556 1.00
## modetectability1[5]     0.04      0.03     0.00     0.11       5494 1.00
## modetectability1[6]     0.02      0.02     0.00     0.07       4917 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).
marginal_effects(m5, categorical = T)

marginal_effects(m5, categorical = F)
## Warning: Predictions are treated as continuous variables in
## 'marginal_effects' by default, which is likely invalid for ordinal
## families. Please set 'categorical' to TRUE.

## Warning: Predictions are treated as continuous variables in
## 'marginal_effects' by default, which is likely invalid for ordinal
## families. Please set 'categorical' to TRUE.

## Warning: Predictions are treated as continuous variables in
## 'marginal_effects' by default, which is likely invalid for ordinal
## families. Please set 'categorical' to TRUE.

Save current workspace to be able to compare which model fits best afterwards

save.image("models.RData")

Analyse which model fits our data best

load("models.RData")

m1 <- add_criterion(m1, "waic")
m2 <- add_criterion(m2, "waic")
m3 <- add_criterion(m3, "waic")
m4 <- add_criterion(m4, "waic")
m5 <- add_criterion(m5, "waic")

loo_compare(m1, m2, m3, m4, m5, criterion = "waic")
##    elpd_diff se_diff
## m4     0.0       0.0
## m2  -143.9      18.1
## m5 -1340.0      49.0
## m3 -1416.9      51.1
## m1 -1417.4      51.0